Analysis of moss change
Wouter Van Landuyt, Hans Van Calster
2021-06-08
## Loading required package: assertthat##
## Attaching package: 'assertthat'## The following object is masked from 'package:tibble':
##
## has_name## Loading required package: glue##
## Attaching package: 'glue'## The following object is masked from 'package:dplyr':
##
## collapse## Loading required package: rlang##
## Attaching package: 'rlang'## The following object is masked from 'package:assertthat':
##
## has_name## The following objects are masked from 'package:purrr':
##
## %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
## flatten_lgl, flatten_raw, invoke, list_along, modify, prepend,
## splice## `summarise()` has grouped output by 'ifbl_4by4', 'Jaar', 'TaxonIDParent'. You can override using the `.groups` argument.##
## -- Column specification --------------------------------------------------------
## cols(
## id = col_double(),
## code = col_character(),
## naam_nederlands = col_character(),
## naam_wetenschappelijk = col_character(),
## taxn_vrs_key = col_character(),
## taxongroep = col_character(),
## usageKey = col_double(),
## scientificName = col_character(),
## rank = col_character(),
## order = col_character(),
## matchType = col_character(),
## phylum = col_character(),
## kingdom = col_character(),
## genus = col_character(),
## class = col_character(),
## confidence = col_double(),
## synonym = col_logical(),
## status = col_character(),
## family = col_character()
## )## i Using '\',\'' as decimal and '\'.\'' as grouping mark. Use `read_delim()` for more control.##
## -- Column specification --------------------------------------------------------
## cols(
## substraat = col_character(),
## substrate = col_character()
## )1 Exploratory data analysis
1.1 General
- Anthocerotophyta = hornworth = hauwmossen
- Bryophyta = bryophytes = bladmossen
- Marchantiophyta = liverworths = levermossen
| phylum | n_obs | n_soorten | n_hokken |
|---|---|---|---|
| Anthocerotophyta | 218 | 4 | 123 |
| Bryophyta | 95624 | 410 | 796 |
| Marchantiophyta | 17486 | 109 | 727 |
## `summarise()` has grouped output by 'phylum'. You can override using the `.groups` argument.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Joining, by = "substraat"1.2 Ellenberg
| phylum | ell_l | number of species |
|---|---|---|
| Bryophyta | 1 | 3 |
| Bryophyta | 2 | 7 |
| Bryophyta | 3 | 14 |
| Bryophyta | 4 | 32 |
| Bryophyta | 5 | 42 |
| Bryophyta | 6 | 78 |
| Bryophyta | 7 | 83 |
| Bryophyta | 8 | 65 |
| Bryophyta | 9 | 8 |
| Marchantiophyta | 1 | 1 |
| Marchantiophyta | 3 | 6 |
| Marchantiophyta | 4 | 11 |
| Marchantiophyta | 5 | 15 |
| Marchantiophyta | 6 | 28 |
| Marchantiophyta | 7 | 25 |
| Marchantiophyta | 8 | 10 |
| phylum | ell_r | number of species |
|---|---|---|
| Bryophyta | 1 | 7 |
| Bryophyta | 2 | 25 |
| Bryophyta | 3 | 40 |
| Bryophyta | 4 | 28 |
| Bryophyta | 5 | 50 |
| Bryophyta | 6 | 64 |
| Bryophyta | 7 | 77 |
| Bryophyta | 8 | 33 |
| Bryophyta | 9 | 8 |
| Marchantiophyta | 1 | 10 |
| Marchantiophyta | 2 | 18 |
| Marchantiophyta | 3 | 9 |
| Marchantiophyta | 4 | 15 |
| Marchantiophyta | 5 | 22 |
| Marchantiophyta | 6 | 13 |
| Marchantiophyta | 7 | 6 |
| Marchantiophyta | 8 | 2 |
| Marchantiophyta | 9 | 1 |
| phylum | ell_f | number of species |
|---|---|---|
| Bryophyta | 1 | 6 |
| Bryophyta | 2 | 6 |
| Bryophyta | 3 | 21 |
| Bryophyta | 4 | 60 |
| Bryophyta | 5 | 73 |
| Bryophyta | 6 | 65 |
| Bryophyta | 7 | 25 |
| Bryophyta | 8 | 37 |
| Bryophyta | 9 | 28 |
| Bryophyta | 10 | 8 |
| Bryophyta | 11 | 1 |
| Bryophyta | 12 | 2 |
| Marchantiophyta | 4 | 8 |
| Marchantiophyta | 5 | 20 |
| Marchantiophyta | 6 | 17 |
| Marchantiophyta | 7 | 18 |
| Marchantiophyta | 8 | 17 |
| Marchantiophyta | 9 | 13 |
| Marchantiophyta | 10 | 2 |
| Marchantiophyta | 11 | 1 |
| phylum | ell_n | number of species |
|---|---|---|
| Bryophyta | 1 | 19 |
| Bryophyta | 2 | 77 |
| Bryophyta | 3 | 63 |
| Bryophyta | 4 | 80 |
| Bryophyta | 5 | 56 |
| Bryophyta | 6 | 28 |
| Bryophyta | 7 | 9 |
| Marchantiophyta | 1 | 19 |
| Marchantiophyta | 2 | 30 |
| Marchantiophyta | 3 | 20 |
| Marchantiophyta | 4 | 14 |
| Marchantiophyta | 5 | 8 |
| Marchantiophyta | 6 | 4 |
| Marchantiophyta | 7 | 1 |
| phylum | ell_t | number of species |
|---|---|---|
| Bryophyta | 2 | 31 |
| Bryophyta | 3 | 117 |
| Bryophyta | 4 | 53 |
| Bryophyta | 5 | 62 |
| Bryophyta | 6 | 42 |
| Bryophyta | 7 | 20 |
| Bryophyta | 8 | 7 |
| Marchantiophyta | 2 | 13 |
| Marchantiophyta | 3 | 37 |
| Marchantiophyta | 4 | 16 |
| Marchantiophyta | 5 | 16 |
| Marchantiophyta | 6 | 9 |
| Marchantiophyta | 7 | 2 |
| Marchantiophyta | 8 | 3 |
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2## Warning: Removed 27 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 27 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 27 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 28 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 58 rows containing missing values## Warning: Removed 27 rows containing missing values (geom_point).## Warning: Removed 16 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 16 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 17 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 47 rows containing missing values## Warning: Removed 27 rows containing missing values (geom_point).## Warning: Removed 16 rows containing missing values (geom_point).## Warning: Removed 16 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 17 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 47 rows containing missing values## Warning: Removed 28 rows containing missing values (geom_point).## Warning: Removed 17 rows containing missing values (geom_point).
## Warning: Removed 17 rows containing missing values (geom_point).## Warning: Removed 17 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 47 rows containing missing values## Warning: Removed 58 rows containing missing values (geom_point).## Warning: Removed 47 rows containing missing values (geom_point).
## Warning: Removed 47 rows containing missing values (geom_point).
## Warning: Removed 47 rows containing missing values (geom_point).## Warning: Removed 47 rows containing non-finite values (stat_density).
## Warning: Removed 4 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 4 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 4 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 4 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 10 rows containing missing values## Warning: Removed 4 rows containing missing values (geom_point).## Warning: Removed 2 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 2 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 2 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 8 rows containing missing values## Warning: Removed 4 rows containing missing values (geom_point).## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 2 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 2 rows containing missing values## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 8 rows containing missing values## Warning: Removed 4 rows containing missing values (geom_point).## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 2 rows containing non-finite values (stat_density).## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 8 rows containing missing values## Warning: Removed 10 rows containing missing values (geom_point).## Warning: Removed 8 rows containing missing values (geom_point).
## Warning: Removed 8 rows containing missing values (geom_point).
## Warning: Removed 8 rows containing missing values (geom_point).## Warning: Removed 8 rows containing non-finite values (stat_density).
1.3 Substrate
Only remove substrate levels with less than 5 species (for Ellenberg variables this is not necessary because they are included as continuous covariate).
| phylum | substrate | number of species |
|---|---|---|
| Bryophyta | artificial stone | 90 |
| Bryophyta | decaying vegetation | 21 |
| Bryophyta | decorticated wood | 48 |
| Bryophyta | epiphytic on living wood | 75 |
| Bryophyta | gravel or sand | 109 |
| Bryophyta | mineral soil | 214 |
| Bryophyta | peat | 48 |
| Bryophyta | rock, hard | 103 |
| Bryophyta | rock, soft | 40 |
| Bryophyta | soil on rock | 69 |
| Marchantiophyta | artificial stone | 10 |
| Marchantiophyta | bryophyte | 20 |
| Marchantiophyta | decaying vegetation | 19 |
| Marchantiophyta | decorticated wood | 16 |
| Marchantiophyta | epiphytic on living wood | 15 |
| Marchantiophyta | gravel or sand | 37 |
| Marchantiophyta | mineral soil | 74 |
| Marchantiophyta | peat | 33 |
| Marchantiophyta | rock, hard | 17 |
| Marchantiophyta | rock, soft | 16 |
| Marchantiophyta | soil on rock | 25 |
| phylum | n | Number of species associated with n substrates |
|---|---|---|
| Bryophyta | 1 | 104 |
| Bryophyta | 2 | 103 |
| Bryophyta | 3 | 61 |
| Bryophyta | 4 | 27 |
| Bryophyta | 5 | 19 |
| Bryophyta | 6 | 9 |
| Bryophyta | 7 | 6 |
| Bryophyta | 8 | 2 |
| Bryophyta | 9 | 1 |
| Marchantiophyta | 1 | 26 |
| Marchantiophyta | 2 | 24 |
| Marchantiophyta | 3 | 12 |
| Marchantiophyta | 4 | 10 |
| Marchantiophyta | 5 | 11 |
| Marchantiophyta | 6 | 9 |
| Marchantiophyta | 7 | 1 |
| Marchantiophyta | 8 | 2 |
1.4 Persisting, lost, gained species
| phylum | trajectory | Number of species |
|---|---|---|
| Bryophyta | gained | 17 |
| Bryophyta | lost | 11 |
| Bryophyta | persisting, decreased | 106 |
| Bryophyta | persisting, increased | 178 |
| Bryophyta | persisting, stable | 20 |
| Marchantiophyta | gained | 2 |
| Marchantiophyta | lost | 12 |
| Marchantiophyta | persisting, decreased | 49 |
| Marchantiophyta | persisting, increased | 28 |
| Marchantiophyta | persisting, stable | 4 |
1.5 Period 1 vs period 2
## Warning: Ignoring unknown aesthetics: text## `geom_smooth()` using formula 'y ~ x'## Warning: Ignoring unknown aesthetics: text## `geom_smooth()` using formula 'y ~ x'| phylum | name | value_binned | n |
|---|---|---|---|
| Bryophyta | n_hok_1980_1999 | [0,3] | 81 |
| Bryophyta | n_hok_1980_1999 | (3,10] | 57 |
| Bryophyta | n_hok_1980_1999 | (10,29] | 64 |
| Bryophyta | n_hok_1980_1999 | (29,83.4] | 65 |
| Bryophyta | n_hok_1980_1999 | (83.4,226] | 65 |
| Bryophyta | n_hok_2000_2019 | [0,3] | 74 |
| Bryophyta | n_hok_2000_2019 | (3,10] | 57 |
| Bryophyta | n_hok_2000_2019 | (10,29] | 57 |
| Bryophyta | n_hok_2000_2019 | (29,83.4] | 63 |
| Bryophyta | n_hok_2000_2019 | (83.4,226] | 81 |
| Marchantiophyta | n_hok_1980_1999 | [0,3] | 19 |
| Marchantiophyta | n_hok_1980_1999 | (3,10] | 19 |
| Marchantiophyta | n_hok_1980_1999 | (10,29] | 20 |
| Marchantiophyta | n_hok_1980_1999 | (29,83.4] | 25 |
| Marchantiophyta | n_hok_1980_1999 | (83.4,226] | 12 |
| Marchantiophyta | n_hok_2000_2019 | [0,3] | 27 |
| Marchantiophyta | n_hok_2000_2019 | (3,10] | 17 |
| Marchantiophyta | n_hok_2000_2019 | (10,29] | 21 |
| Marchantiophyta | n_hok_2000_2019 | (29,83.4] | 17 |
| Marchantiophyta | n_hok_2000_2019 | (83.4,226] | 13 |
2 Methods
2.1 Trends
To assess the changes in distribution of bryophyte species form 1980 until 2019 we divided the survey in two periods and compared the period 1980-1999 with the period 2000-2019. To minimize the effect of badly prospected grid cells we only used grid cells (16km²) where at least 40 species were found in both periods (Figure 2.1).
Figure 2.1: 16 km² grids of Flanders and the Brussels Capital Region with at least 40 species of Bryophytes recorded in both the period 1980-1999 and the period 2000-2019.
This resulted in 227 grid cells prospected in both periods. In these selected grid cells 35648 bryophyte observations (expressed as unique species per year and per 16 km² grid) were done in the first period and 45236 in the second period. To correct for potential changes in survey effort within those grids the number of occupied grid cells in the second period was modeled as a function of the number of occupied grid cells in the first period so that effects of other explanatory variables could be interpreted relative to the change in survey effort. The conceptual idea comes from Telfer, Preston, and Rothery (2002), but the technical implementation that we followed is different from the original proposal.
The method assumes that during the two surveys recorders attempted to observe as many species as possible in each grid cell. The surveyors in both periods looked for liverworths and mosses and we therefor did not consider survey effort to differ between both phyla. The size of the distribution area of individual species in a given period, i.e., the number of 16 km² grid cells where the species was observed, was expressed as a proportion of the total number of grid cells studied during that period. To summarise the changes in distribution area for all species, we modeled the number of occupied 16 km² squares in the second survey (out of the total number of 16 km² squares considered, i.e. proportional data) as a function of the logit-transformed proportion in the first survey (baseline) and covariates (phylum, substrate and Ellenberg values for light, moisture, nitrogen and temperature). Substrate effects (11 levels) were modeled for each phylum separately as random effects coming from a normal distribution with mean zero and an estimated standard deviation. A species could be associated with more than substrate type. If that was the case, a weighting factor (one divided by the number of substrate types) was used to ensure that each species in a given period has equal weight in the analysis (if this is not done, confidence interval would be too narrow and estimates could be biased). Because the counts were overdispersed (variance larger than Binomial distribution), a beta-binomial distribution was used to account for overdispersion and a logit-link was used to link the proportions to the linear predictor with the following model equation:
\[ g(\mu_{i,s2}) = \beta_0 + \beta_1\log\frac{\mu_{i,s1}}{1-\mu_{i,s1}} + (\beta_L\text{ell}_L + \beta_F\text{ell}_F + \beta_N\text{ell}_N + \beta_T\text{ell}_T)\text{phylum}_k + \mathbf{b}_{o,j} + \mathbf{b}_{1,j}\text{phylum}_k \] where:
\[ \begin{aligned} \mu_{i,s2} \equiv \mathbf{E}(Y_{i,s2}|\text{trials, weights},\mathbf{b})&& \text{expected value}\\ g(\mu_{i,s2}) = \log\frac{\mu_{i,s2}}{1-\mu_{i,s2}} && \text{logit-link} \\ Y_{i,s2} \sim \text{BetaBinomial}(\mu_{i,s2}, \phi_{i,s2}) && \text{Beta-Binomial distribution}\\ \begin{bmatrix}\mathbf{b}_{o,j}\\ \mathbf{b}_{1,j}\end{bmatrix}\sim N(0,\Sigma) && \text{random intercept and slope}\\ \Sigma = \begin{bmatrix}\sigma^2_{0,j}&0\\0&\sigma^2_{1,j}\end{bmatrix} && \text{variance-covariance matrix} \end{aligned} \]
All species were included also those that occurred in less than five 16 km² grid cells in the first period. In Telfer, Preston, and Rothery (2002), these species were removed, but our model diagnostics showed that they could be included. To avoid calculating the log of zero, the proportion of occupied grid cells in the first period, \(\mu_{i,s1}\) was calculated as \((n_{i, s1}+0.5)/(n_{\text{tot}}+1)\). This transformation is not necessary for the response variable because the beta-binomial accommodates zeroes.
A Bayesian hierarchical modeling procedure was used to calculate the posterior distribution for these parameters (Bürkner 2017) and we assumed non-informative priors.
We focused model selection on the Ellenberg terms only. All other effects were kept in the model. First, interaction terms for Ellenberg-values with phylum for which the 95% confidence interval contained 0 were removed. Next, main Ellenberg terms for which the 95% confidence interval contained 0 were removed. The simplest model with the highest predictive accuracy based on leave-one-out cross validation (Vehtari, Gelman, and Gabry 2017) was selected.
The Ellenberg terms are in principle ordinal values. In an exploratory phase, we modeled these terms as monotonic effects (Bürkner and Charpentier 2020), but the results reassured us that we could treat the Ellenberg values as continuous values.
All analyses were done using the R language for statistical computing (R Core Team 2020).
3 Results
3.1 Model selection
| elpd_diff | se_diff | elpd_loo | se_elpd_loo | p_loo | se_p_loo | looic | se_looic | |
|---|---|---|---|---|---|---|---|---|
| alt_telfer_2b | 0.00 | 0.00 | -4178.25 | 48.38 | 62.88 | 8.36 | 8356.50 | 96.77 |
| alt_telfer_2a | -3.96 | 1.75 | -4182.21 | 48.71 | 69.58 | 9.24 | 8364.42 | 97.42 |
| alt_telfer | -10.22 | 2.79 | -4188.48 | 48.50 | 78.43 | 9.69 | 8376.95 | 97.00 |
3.2 Model diagnostics
plot(loo1$loos$alt_telfer_2b, label_points = TRUE)
## Using 10 posterior samples for ppc type 'dens_overlay' by default.
3.3 Model summary of best predictive model
## Warning in tidy.brmsfit(alt_telfer_2b): some parameter names contain
## underscores: term naming may be unreliable!| effect | component | group | term | estimate | std.error | conf.low | conf.high |
|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 0.5041175 | 0.1786289 | 0.1524429 | 0.8514066 |
| fixed | cond | NA | qlogisn_hok_1980_1999P0.5Dn_hok_tot_1980_1999P1 | 0.8606107 | 0.0244643 | 0.8119928 | 0.9097005 |
| fixed | cond | NA | ell_n | 0.0594435 | 0.0306046 | -0.0005739 | 0.1190961 |
| fixed | cond | NA | ell_f | -0.1036722 | 0.0228619 | -0.1476640 | -0.0584242 |
| fixed | cond | NA | ell_n:phylum1 | -0.0733091 | 0.0276106 | -0.1268828 | -0.0184745 |
| fixed | cond | NA | ell_f:phylum1 | 0.0619373 | 0.0157723 | 0.0310884 | 0.0931486 |
| ran_pars | cond | substrate | sd__(Intercept) | 0.3000726 | 0.0970315 | 0.1661444 | 0.5352504 |
| ran_pars | cond | substrate | sd__phylum1 | 0.0746576 | 0.0607954 | 0.0029444 | 0.2277854 |
3.4 Model visualisation
3.4.1 Y-scale = expected number of occurrences (out of 228 grid cells)
Figure 3.1:
Substrate random effects are calculated from the posterior linear predictor, conditional on Ellenberg values and the number of occupied grid cells fixed at their mean values. Gelman, Hill, and Yajima (2012) point out that no correction for multiple comparisons are necessary when substrate effects are calculated in this way.
Figure 3.1: Estimated total number of occupied cells. The total number of grid cells in both periods equals 228. Expectations and 95% confidence bounds or intervals are given. These expectations are conditional on n_hok_1980_1999 = 52.04 & ell_f = 5.87 & ell_n = 3.4 for variables other than the one on the x-axis (held constant to mean value). The black point with 95% confidence ranges represents the overall effect conditional on all other variables and can be used as a reference across all figures. (a) Change in recording probability. Black line: line of no change. (b) Effect of substrate on occurrence of mosses. (c) Effect of Ellenberg nitrogen values on occurrence of mosses. (d) Effect of Ellenberg moisture values on occurrence of mosses.